** This .do file can be run separately, but feeds into the main DHSV (bootstrap) file to produce the bootstrapped standard errors reported in the paper.

** Set file paths to folders to store intermediate data files, output and .do files:
global datafiles "N:\SpecialProjects\Inequality_Egypt\Peru_2019\data"
global dofiles "N:\SpecialProjects\Inequality_Egypt\Peru_2019\do_files"

** Policy and analysis years for main specification (as seen in Table 4 (and B.3 and B.4))
local aby "1990"
local aey "1997"
local pby "1996"
local pey "1997"
local exc "& year~=1995"

/*
** alternative specifications (as seen in Table 6):

** extend policy years from 1995 to 1998:
local aby "1990"
local aey "1998"
local pby "1995"
local pey "1998"
local exc " "

** extend policy years from 1995 to 1999:
local aby "1990"
local aey "1999"
local pby "1995"
local pey "1999"
local exc " "

*/



************* Logit (as described in Section 6 and shown in Appendix Table B.1) ******************************

use "$datafiles/peru_dhsV_analysis.dta", clear

keep caseid age_* ster_* kids_*  mar_* will_end*  boys_*  born_* m_age_fb yr_st ed_cat weight1 rur_coast urb_mount rur_mount urb_jungle rur_jungle urb_coast  ///

*Reshape to create a pseudo-panel
reshape long age_ ster_ kids_  mar_ will_end_  boys_  born_  , i(caseid) j(year)


* dummy for being sterilized in a given year
gen sterilized=yr_st==year


* Drop observations after sterilization (no longer at risk of sterilization)
bys caseid: gen flag=cond(sterilized==0, sum(sterilized),0)
drop if flag==1


* It was "illegal" to be sterilized if you were not married and did not have kids - and almost no one was, so we make this 
* sample restriction (otherwise the probit does not converge due to perfect prediction of failure.)
gen eligable=(mar_ | will_end_) & kids_>0


gen kids_count=kids_
replace kids_count=. if kids_==0
replace kids_count=10 if kids_>10
tab kids_count, gen(kids)

gen year_lin=year-1992

gen young=age_<26
gen mid1_age=age_>=26 & age_<=29
gen mid2_age=age_>29 & age_<=36
gen older=age_>36

gen ed_prim=ed_cat==2 | ed_cat==1
gen ed_second=ed_cat==3
gen ed_high=ed_cat==4

gen no_birth=born_==0

gen no_boys=boys_==0


local int_vars " rur_mount rur_jungle ed_prim young mid1_age no_birth "


*** (Appendix Table B.1)
local r "replace"
foreach p of numlist 1 2  {

** Actual
local P1  "year==1996 | year==1997"
local Y1  "year>=1990 & year<=1997 & year~=1995"
** Placebo
local P2  "year==1993 | year==1994"
local Y2  "year>=1990 & year<=1994" 

cap drop policy 
gen policy=0
replace policy=1 if `P`p''


cap drop policy_X_*
foreach x of local int_vars{
qui gen policy_X_`x'=`x'*policy
}


forval i=1/9{
gen policy_X_kids`i'=policy*kids`i'
}


	logit sterilized  kids1-kids9 m_age_fb      rur_coast  urb_mount     rur_mount urb_jungle rur_jungle ed_prim ed_second young mid1_age older  no_boys no_birth ///
       policy_X_rur_mount policy_X_rur_jungle policy_X_ed_prim policy_X_young policy_X_mid1_age policy_X_no_birth  year_lin [pw=weight1] if `Y`p'' & eligable, vce(cluster caseid)
	
	outreg2 using "$datafiles/logit_V_placebos",  excel  bdec(3) `r' ctitle("P =`P`p''","Y=`Y`p''") 
	
	local r "append"
}	


local int_vars " rur_mount rur_jungle ed_prim young mid1_age no_birth "

cap drop policy
gen policy=0
replace policy=1 if year>=`pby' & year<=`pey'


cap drop policy_X_*
foreach x of local int_vars{
qui gen policy_X_`x'=`x'*policy
}

forval i=1/9{
gen policy_X_kids`i'=policy*kids`i'
}

cap drop policy_save
gen policy_save=.
foreach x of local int_vars{
qui gen policy_X_`x'_save=.
}

forval i=1/9{
gen policy_X_kids`i'_save=.
}
	

	
logit sterilized  kids1-kids9 m_age_fb rur_coast urb_mount rur_mount urb_jungle rur_jungle ed_prim ed_second young mid1_age older  no_boys no_birth  ///
      policy_X_rur_mount policy_X_rur_jungle policy_X_ed_prim policy_X_young policy_X_mid1_age  policy_X_no_birth  year_lin [pw=weight1] if year>=`aby' & year<=`aey' `exc' & eligable, vce(cluster caseid)


qui cap drop pr2_

qui predict pr2_ if year>=`pby' & year<=`pey' `exc' & eligable & e(sample)


qui replace policy_save=policy
foreach x of local int_vars{
qui replace policy_X_`x'_save=policy_X_`x'
}

forval i=1/9{
qui replace policy_X_kids`i'_save=policy_X_kids`i'
}


qui replace policy=0
foreach x of local int_vars{
qui replace policy_X_`x'=0
}

forval i=1/9{
qui replace policy_X_kids`i'=0
}

qui cap drop pr1_

qui predict pr1_ if year>=`pby' & year<=`pey' `exc'  & eligable & e(sample)


qui replace policy=policy_save
foreach x of local int_vars{
qui replace policy_X_`x'=policy_X_`x'_save
}

forval i=1/9{
qui replace policy_X_kids`i'=policy_X_kids`i'_save
}

save "$datafiles/peru_dhsV_qp.dta", replace

** All de jure women who were eiligable in 1996 or 1997 have a pr1 pr2 for that year 
** Being eligible means having ever been married and having at least one child
keep pr1_* pr2_* caseid year 

*eligable sterilized

* make this wide again
*reshape wide eligable sterilized pr1_ pr2_  , i(caseid) j(year)

reshape wide pr1_ pr2_  , i(caseid) j(year)
sort caseid
save  "$datafiles/pscoresV.dta", replace

** Merge the p1 p2 variables onto the cross-sectional data set.
use "$datafiles/peru_dhsV_analysis.dta", clear
cap drop _merge
sort caseid
merge caseid using "$datafiles/pscoresV.dta"

cap drop S
gen S = yr_st>=`pby' & yr_st<=`pey' 


gen status=.
replace status=1 if yr_st<`pby' & yr_st~=.     /* sterilized before the policy -- these women are all excluded from the outcome results below */


egen _p2=rowmean(pr2_*)
egen _p1=rowmean(pr1_*)

gen delta_p=_p2-_p1
gen dp_p2=delta_p/_p2

*** Standard Reweighting 

gen sdrweight=weight1 if S==1 & status~=1
replace sdrweight=weight1*(_p2/(1-_p2)) if S==0 & status~=1

bys S: egen sdsrweight_t=total(sdrweight)
gen sdsrweight=sdrweight/sdsrweight_t


*** Our reweighting S2=1, D=1
gen rweight=weight1*dp_p2 if S==1 & status~=1                 
replace rweight=weight1*(delta_p/(1-_p2)) if S==0 & status~=1 

bys S: egen srweight_t=total(rweight)
gen srweight=rweight/srweight_t

*** Our reweighting S2=1, D=0
cap drop orweight osrweight_t osrweight

gen orweight=weight1*(_p1/_p2) if S==1 & status~=1          
replace orweight=weight1*((_p1)/(1-_p2)) if S==0 & status~=1 

bys S: egen osrweight_t=total(orweight)
gen osrweight=orweight/osrweight_t


save "$datafiles/dataV_w_weights.dta", replace


**************************** Estimating Results *********************************************************************************

*********** Estimating the number of fewer children born due to the campaign ******************************

* In order to estimate the number of women affected we need to scale the weights by the (approximate) number of women each survey respondent represents in the relevant surveyed Peruvian population (women 15-49)
* 325.8 = number of women 15-49 (in the relevant year from UN population tables) divided by the full sample size in DHS
gen w3=weight1*325.8

total w3 if S==1
matrix total=e(b)
scalar women_all=total[1,1]

gen estim=w3*dp_p2
total estim if S==1
matrix total=e(b)
scalar women_campaign=total[1,1]
local women=total[1,1]


**************************** Number of Kids *****************************************************************
* Modified Reweighting (S2=1, C=1) 
reg kids_2004 S  [pw=srweight] if status~=1 
scalar kids_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, replace  excel ctitle("kids", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'"  ) bdec(3) adds(Num Women,`women')  

local r `append'

* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
reg kids_2004 S [pw=osrweight] if status~=1 
scalar kids_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel `r' ctitle("kids", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)


**************************** Mother's Health *****************************************************************

              ** weight_height_m


* Modified Reweighting (S2=1, C=1) 
reg weight_height_m S  [pw=srweight] if status~=1 & pregnant==0
scalar WH_M_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, excel  `r'  ctitle("WH_np", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'") bdec(3)

* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
reg weight_height_m S [pw=osrweight] if status~=1  & pregnant==0
scalar  WH_M_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel  `r'  ctitle("WH_np", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'") bdec(3)


              ** height_age_m


* Modified Reweighting (S2=1, C=1) 
reg height_age_m S  [pw=srweight] if status~=1 & pregnant==0
scalar HA_M_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, excel  `r'  ctitle("HA_np", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)

* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
reg height_age_m S [pw=osrweight] if status~=1  & pregnant==0
scalar HA_M_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel  `r'  ctitle("HA_np", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)


               **Anemic

* Modified Reweighting (S2=1, C=1) 
reg anemic_m S  [pw=srweight] if status~=1 & pregnant==0
scalar anemic_M_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, excel  `r'  ctitle("anemic_np", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'"  ) bdec(3)

* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
reg anemic_m S [pw=osrweight] if status~=1 & pregnant==0
scalar anemic_M_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel  `r'  ctitle("anemic_np", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'"  ) bdec(3)




*********** Domestic Violence LAST 12 MONTHS*******************************************************

* Modified Reweighting (S2=1, C=1) 
reg p_physexviol12m S [pw=srweight] if status~=1
scalar dm_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, excel append ctitle("dm", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'"  ) bdec(3)

* Modified Reweighting (S2=1, C=0) 
reg p_physexviol12m S [pw=osrweight] if status~=1 
scalar dm_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel append ctitle("dm", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'"  ) bdec(3)




**************************** Mother's Labor-Force Participation (worked for pay) *************************************

* Modified Reweighting (S2=1, C=1) 
reg working_p S  [pw=srweight] if status~=1 
scalar wp_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, excel  `r'  ctitle("work_p","`aby' to `aey'", "policy `pby' to `pey'", "`exc'"  ) bdec(3)

* Modified Reweighting (S2=1, C=0) 
reg working_p S [pw=osrweight] if status~=1 
scalar wp_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel  `r'  ctitle("work_p", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)





***********CHILDREN's EDUCATION

**********Education of household children (under 15) but born prior to the program
** need to use the household record to get these kids' education which leads to slight differences


use "$datafiles/dataV_w_weights.dta", clear


*drop kby_*
*drop if caseid==""

keep kid_age* kby* kid_boy_ed* enrolled* sch_yr_* grade_at_age* caseid status S srweight osrweight

reshape long kid_age_ed kbyed kid_boy_ed enrolled sch_yr_ grade_at_age, i(caseid) j(kid)

rename kid_boy_ed kid_boy
rename kid_age_ed kid_age
rename kbyed kby


*********** Education of own children**********************************************************************


*********** Enrollment of own kids**********************************************************************

* Modified Reweighting (S2=1, C=1) 
reg enrolled S i.kid_age  [pw=srweight]        if status~=1 & kby<=`pey'  & kid_age<=15 & kid_age>=6 
scalar ENkid_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, excel append ctitle("enrolled_gt6",  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'"  ) bdec(3)




* Modified Reweighting (S2=1, C=0) 
reg enrolled S i.kid_age [pw=osrweight]        if status~=1 & kby<=`pey'  & kid_age<=15 & kid_age>=6 
scalar ENkid_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel append ctitle("enrolled_gt6", "`aby' to `aey'", "policy `pby' to `pey'", "`exc'"  ) bdec(3)

gen ed_flag=e(sample)

*********** grade for age of own children**********************************************************************

* Modified Reweighting (S2=1, C=1) 
qui reg grade_at_age S i.kid_age  [pw=srweight] if status~=1 & kby<=`pey' & kid_age<=15 & kid_age>=6 & ed_flag==1
scalar GFAkid_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, excel `r' ctitle("grade_f_age",  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'"  ) bdec(3)

* Modified Reweighting (S2=1, C=0) 
qui reg grade_at_age S i.kid_age  [pw=osrweight] if status~=1 & kby<=`pey'  & kid_age<=15 & kid_age>=6  & ed_flag==1 
scalar GFAkid_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel `r' ctitle("grade_f_age",  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'"  ) bdec(3)


**********************ED by GENDER ***************************************************************************


*********** Enrollment of own boys**********************************************************************

* Modified Reweighting (S2=1, C=1) - only positive weigths are included by default since pw drops negative weights
reg enrolled S i.kid_age  [pw=srweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>4 & kid_boy==1
scalar ENboy_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, excel append ctitle(boy_enrolled ,  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)


* Modified Reweighting (S2=1, C=0) - only positive weigths are included by default since pw drops negative weights
reg enrolled S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>4 & kid_boy==1 
scalar ENboy_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel append ctitle(boy_enrolled,  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3) 


*********** Enrollment of own girls**********************************************************************

* Modified Reweighting (S2=1, C=1) 
reg enrolled S i.kid_age  [pw=srweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>4 & kid_boy==0
scalar ENgirl_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, excel append ctitle(girl_enrolled,  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)


* Modified Reweighting (S2=1, C=0) 
reg enrolled S i.kid_age [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>4 & kid_boy==0 
scalar ENgirl_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel append ctitle(girl_enrolled,  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)


 
********** grade for age of own boys **********************************************************************

* Modified Reweighting (S2=1, C=1) 
qui reg grade_at_age S i.kid_age  [pw=srweight] if status~=1 & kby<=1997 & kid_age<=15 & kid_age>=6  & kid_boy==1  & ed_flag==1
scalar GFAboy_s2d1=_b[S]
outreg2 using $datafiles\DHS_V_results, excel `r' ctitle(grade_f_age,  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)

* Modified Reweighting (S2=1, C=0) 
qui reg grade_at_age S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6 & kid_boy==1  & ed_flag==1
scalar GFAboy_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel `r' ctitle(grade_f_age,  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)

********** grade for age of own girls **********************************************************************

* Modified Reweighting (S2=1, C=1) 
qui reg grade_at_age S i.kid_age  [pw=srweight] if status~=1 & kby<=1997 & kid_age<=15 & kid_age>=6  & kid_boy==0  & ed_flag==1
scalar GFAgirl_s2d1=_b[S] 
outreg2 using $datafiles\DHS_V_results, excel `r' ctitle(grade_f_age,  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)

* Modified Reweighting (S2=1, C=0) 
qui reg grade_at_age S i.kid_age  [pw=osrweight] if status~=1 & kby<=1997  & kid_age<=15 & kid_age>=6  & kid_boy==0  & ed_flag==1
scalar GFAgirl_s2d0=_b[S]
outreg2 using $datafiles\DHS_V_results, excel `r' ctitle(grade_f_age,  "`aby' to `aey'", "policy `pby' to `pey'", "`exc'" ) bdec(3)


matrix non_boot=(kids_s2d1,	kids_s2d0, WH_M_s2d1, WH_M_s2d0, HA_M_s2d1, HA_M_s2d0, anemic_M_s2d1, anemic_M_s2d0,  wp_s2d1,	 wp_s2d0, dm_s2d1, dm_s2d0,	 ENkid_s2d1, ENkid_s2d0 , GFAkid_s2d1, GFAkid_s2d0, ENboy_s2d1, ENboy_s2d0,  ENgirl_s2d1, ENgirl_s2d0, GFAboy_s2d1, GFAboy_s2d1, GFAgirl_s2d1, GFAgirl_s2d1, women_all, women_campaign )
 
matrix list non_boot



************ blancing table (Appenddix Table B.2) **************************************************
use "$datafiles/dataV_w_weights.dta", clear

gen var_1=age_1996
gen var_2=kids_1996
gen var_3=ed_years
gen var_4=m_age_fb
gen var_5=rural
gen var_6=sreg1
gen var_7=sreg2
gen var_8=sreg3
gen var_9=rtype1
gen var_10=rtype2
gen var_11=rtype3
gen var_12=rtype4
gen var_13=rtype5
gen var_14=rtype6


cap drop flag1
reg num_kid S [pw=srweight] if status~=1
gen flag1=0
replace flag1=1 if e(sample)


* make a file to hold the table
preserve
	clear
	set obs 14
	gen colvar=""
	forval v=1/14 {
	replace colvar="var_`v'" in `v'
				}
				forval j=1/5{
		gen specification_`j'=.
	}
gen pval=.		

	save "$datafiles/datatable_balanceDHSV.dta", replace emptyok
restore

		

* pvalue of differences under our preferred reweighting
forvalues v=1/14 {
	qui svyset [pweight=srweight]
	qui svy: mean var_`v' , over(S)
	qui lincomest [var_`v']0 - [var_`v']1
	matrix results=r(table)
	local pvalue_`v'=results[4,1]	
	preserve
		use "$datafiles/datatable_balanceDHSV.dta", clear
		qui replace pval=`pvalue_`v''  in  `v'
		save "$datafiles/datatable_balanceDHSV.dta", replace
		restore
	}

local spec_1 "if dp_p2~=. & S==1 [aw=weight1]"    /* not reweighted - sterilized */
local spec_2 "if dp_p2~=. & S==1 [aw=srweight]"   /* reweighted - program, sterilized */
local spec_3 "if dp_p2~=. & S==1 [aw=osrweight]"  /* reweighted - non program, sterilized */
local spec_4 "if dp_p2~=. & S==0 [aw=weight1]"    /* not reweighted - not sterilized */
local spec_5 "if dp_p2~=. & S==0 [aw=srweight]"   /* reweighted - program, not sterilized */


forvalues j=1/5{
	forvalues v=1/14{
		qui sum var_`v' `spec_`j''
		local table_`v'_`j'=r(mean)
		preserve
		use "$datafiles/datatable_balanceDHSV.dta", clear
		qui replace specification_`j'=`table_`v'_`j''  in  `v'
		save "$datafiles/datatable_balanceDHSV.dta", replace
		restore
		}
}



use "$datafiles/datatable_balanceDHSV.dta", clear

outsheet using "$datafiles/Balance_DHSV_2020.xls", replace






